# see \docs\audio\uade 2.13\contrib
# see \docs\audio\docs\UADE - Developer discussion
# see \docs\audio\modguide
# tested with python 3.2, numpy 1.6.1, scipy 0.10.1

import math
import csv;
from numpy import linspace
from scipy import exp, array
from scipy.signal import firwin, lsim
from scipy.fftpack import fft, ifft

AMIGA_PAL_CLOCK = 3546895
VIDEO_MODE = 'pal'
AMIGA_PAL_CLOCK = 3579545 # ntsc
     
def log_safe(x):
    if (x <= 0):
        return float('-1e-100')
    else:
        return math.log(x);
    
def fir_minphase(table):
    
    table = list(table)
    
    # table should be a real-valued table of FIR coefficients
    convolution_size = len(table)
    
    # pad table with zeros
    table += [0] * (convolution_size * 7)
    
    # compute the real cepstrum
    # fft -> abs + ln -> ifft
    cepstrum = ifft([log_safe(x) for x in abs(fft(table))])

    # window the cepstrum in such a way that anticausal components become rejected
    cepstrum[1:len(cepstrum)/2] *= 2;
    cepstrum[len(cepstrum)/2+1:len(cepstrum)  ] *= 0;
    
    # now cancel the previous steps:
    # fft -> exp -> ifft -> real
    cepstrum = ifft([exp(x) for x in fft(cepstrum)])
    
    return [x.real for x in cepstrum[0:convolution_size]]

def quantize(xx):
    
    x = list(xx)
    fact = 2 ** 17
    
    # this adjusts range precisely between -65536 and 0 so that our bleps look right.
    # we should only do this if the table is integrated; in fact this code probably
    # belonged to the integration function.
    correction_factor = x[-1] - x[0]

    for i in range(len(xx)):
        val = x[i] * fact / correction_factor;
        # correct rounding
        if val < 0:
            val = int(val - 0.5)
        else:
            val = int(val + 0.5)
        x[i] = -val;
    return x
            
def write_bleps(tables, filename, header):

    with open(filename, 'w') as f:
        f.write(header);
        f.write('\n');
        for i in range(len(tables[1])):
            f.write("{0},".format(i));
            for t in tables:
                f.write("{0},".format(t[i]))
            f.write("\n");
            
def write_csv(table, filename):

    with open(filename, 'w') as f:
        for i in range(len(table)):
            f.write("{0}\n".format( table[i]))

def step_calc(table):

    startval = -sum(table)
    result = []
    for el in table:
        startval += el;
        result.append(startval)
    
    return result

def prepare_unfiltered():

    # Because Amiga only has 84 dB SNR, the noise floor is low enough with -90 dB.
    # A500 model uses slightly lower-quality kaiser window to obtain slightly
    # steeper stopband attenuation. The fixed filters attenuates the sidelobes by
    # 12 dB, compensating for the worse performance of the kaiser window.

    # 21 kHz stopband is not fully attenuated by 22 kHz. If the sampling frequency
    # is 44.1 kHz, all frequencies above 22 kHz will alias over 20 kHz, thus inaudible.
    # The output should be aliasingless for 48 kHz sampling frequency.
    a500_unfiltered  = firwin(2048, 21000.0 / AMIGA_PAL_CLOCK * 2, window=('kaiser', 8.0))
    a1200_unfiltered = firwin(2048, 21000.0 / AMIGA_PAL_CLOCK * 2, window=('kaiser', 9.0))
    
    write_csv(a500_unfiltered, 'sinc-A500-{0}.csv'.format(VIDEO_MODE))
    write_csv(a1200_unfiltered, 'sinc-A1200-{0}.csv'.format(VIDEO_MODE))
    
    # move filtering effects to start to allow IIRs more time to settle
    a500_unfiltered  = fir_minphase(a500_unfiltered)
    a1200_unfiltered = fir_minphase(a1200_unfiltered)
    
    write_csv(a500_unfiltered, 'sinc-min-phase-A500-{0}.csv'.format(VIDEO_MODE))
    write_csv(a1200_unfiltered, 'sinc-min-phase-A1200-{0}.csv'.format(VIDEO_MODE))

    return a500_unfiltered, a1200_unfiltered
            
def blep_calc(kind, filter_off, filter_on, unfiltered):
    
    # integrate to produce blep
    filter_off_step = step_calc(filter_off)
    filter_on_step = step_calc(filter_on)
    unfiltered_step = step_calc(unfiltered)
    
    write_bleps([filter_off_step, 
                 filter_on_step, 
                 unfiltered_step], 
                'bleps-unquantized-{0}-{1}.csv'.format(kind, VIDEO_MODE), 
                '"filter off","filter on","unfiltered"')
    
    # quantize and scale
    filter_off_quantize = quantize(filter_off_step)
    filter_on_quantize = quantize(filter_on_step)
    unfiltered_quantize = quantize(unfiltered_step)
    
    write_bleps([filter_off_quantize, 
                 filter_on_quantize,  
                 unfiltered_quantize], 
                'bleps-{0}-{1}.csv'.format(kind, VIDEO_MODE), 
                '"filter off","filter on","unfiltered"')